cd "e:\论文相关"
/*
import excel "datadofile\GDP.xlsx", sheet("Sheet1") firstrow clear
drop if _n<=2
destring  expectGDPgrowth,replace
duplicates drop city year,force
sort province city year
save expectGDPgrowth,replace


import excel "newdata\中国地级市增长目标2000-2019.xlsx", sheet("地区实际GDP增长率") firstrow clear
drop if _n<=2
destring _all,replace
duplicates drop city year,force
sort city year
save RealGDPgrowth,replace

import excel "newdata\中国地级市增长目标2000-2019.xlsx", sheet("Sheet1") firstrow clear
drop if _n<=2
keep 地区 城市名称
rename 地区 districtname
rename 城市名称 city
gen east=.
replace east=1  if (district=="东部") & !missing(district)
replace east=0  if (district!="东部") & !missing(district)
destring _all,replace
duplicates drop city ,force
sort city
save district,replace

import excel "newdata\中国省级增长目标2000-2019.xlsx", sheet("2000-2018") firstrow clear
drop if _n<=2
keep 年份 省份 经济增长目标 
rename 年份 year
rename 省份 province
rename 经济增长目标 provinceexpectGDPgrowth
destring _all,replace
duplicates drop province year,force
sort province year
save provincerxpectGDP,replace

import excel "newdata\1997-2017年中国地级市碳排放数据.xlsx", sheet("数据") firstrow clear
drop if _n<=2
keep CityName year CO2
rename CityName city
destring _all,replace
duplicates drop city year,force
sort city year
save CO2,replace

import excel "控制变量搜集的数据\2000~2020年全国各城市夜间灯光数据.xlsx", sheet("Sheet1") firstrow clear
drop if _n<=2
keep 年份 城市 总和 
rename 年份 year
rename 城市 city
rename 总和 citylight
destring _all,replace
duplicates drop city year,force
sort city year
save citylight,replace

import excel "控制变量搜集的数据\分城市人口情况.xlsx", sheet("sheet1") firstrow clear
drop if _n<=2
keep 年度标识 城市名称 年末总人口万人 城镇个体劳动者万人
rename 年度标识 year
rename 城市名称 city
rename 年末总人口万人 population
rename 城镇个体劳动者万人 citypopulation
destring _all,replace
duplicates drop city year,force
sort city year
save population,replace


import excel "控制变量搜集的数据\汇率.xlsx", sheet("sheet1") firstrow clear
drop if _n<=2
keep 汇率日期 人民币元对美元汇率
rename 汇率日期 date
rename 人民币元对美元汇率 exchangerate
gen year=substr(date,1,4)
destring date,replace ignore("-")
destring _all,replace
bysort year: egen long maxdate=max(date)
keep if date==maxdate
keep year exchangerate
sort year
save exchange,replace



import excel "控制变量搜集的数据\分城市外商投资.xlsx", sheet("sheet1") firstrow clear
drop if _n<=2
keep 年度标识 城市名称 外商实际投资额万美元
rename 年度标识 year
rename 城市名称 city
rename 外商实际投资额万美元 FDI
destring _all,replace
sort  year
merge m:1 year using exchange
duplicates drop city year,force
replace FDI=FDI*exchangerate/10^4
save FDI,replace

import excel "控制变量搜集的数据\分城市国内生产总值.xlsx", sheet("sheet1") firstrow clear
drop if _n<=2
keep 年度标识 城市名称 国内生产总值亿元 第二产业占GDP比重 第三产业占GDP比重
rename 年度标识 year
rename 城市名称 city
rename 国内生产总值亿元 GDP
rename 第二产业占GDP比重 industry2GDP
rename 第三产业占GDP比重 industry3GDP
destring _all,replace
duplicates drop city year,force
sort city year
save industry,replace


import excel "控制变量搜集的数据\分城市固定资产投资.xlsx", sheet("sheet1") firstrow clear
drop if _n<=2
keep 年度标识 城市名称 固定资产投资总额万元
rename 年度标识 year
rename 城市名称 city
rename 固定资产投资总额万元 fixasset
replace fixasset=fixasset/10^4
destring _all,replace
duplicates drop city year,force
sort city year
save fixasset,replace


import excel "newdata\数据求助\市政府工作报告词频分析.xlsx", sheet("2003-2019年各省政府工作报告环境词汇词频分析") firstrow clear
rename 城市 city 
rename 年份 year 
bysort city year: egen numwords=total(词频)
duplicates drop city year,force
keep city year numwords
destring year,replace
sort city year
save numwords,replace

use windspeed
rename 省 province
rename 市 city
rename 年份 year
duplicates drop city year,force
save windspeed,replace

use 20211210
merge 1:1 city year using windspeed,nogen
keep if year>=2000  & year<=2017
save 20211210, replace



clear
import excel "各省地级市数量.xlsx", sheet("Sheet1") firstrow
sort province
save cityquantity,replace
*/
cd "e:\论文相关"
use expectGDPgrowth,clear
duplicates drop city year,force
merge 1:1 city year using CO2,nogen
merge 1:1 city year using RealGDPgrowth,nogen
merge 1:1 city year using citylight,nogen
merge 1:1 city year using industry,nogen
merge 1:1 city year using FDI,nogen
merge 1:1 city year using fixasset,nogen
merge 1:1 city year using population,nogen
merge 1:1 city year using windspeed,nogen
cap drop _m
sort city 
merge     city      using district 
merge m:1 province year using provincerxpectGDP,nogen
/*
cap drop _m
sort province yea
merge m:1 province year using provincerxpectGDP,nogen
*/
keep if year>=2000  & year<=2017
destring RealGDPgrowth,replace force
destring expectGDPgrowth,replace force
destring citylight,replace force
destring GDP,replace force
gen realGDPgrowth=RealGDPgrowth/100
gen CO2perGDP=CO2/GDP
gen FDIGDP=FDI/GDP
replace fixasset=fixasset/GDP

egen Cityname=group(city)

duplicates drop city year,force
xtset Cityname year
*xtdes
gen pressure = ln(expectGDPgrowth/L.realGDPgrowth)
gen pressure2 = expectGDPgrowth/L.realGDPgrowth
gen citylightgrowth= citylight/L.citylight-1
gen lastrealGDP=ln(L.GDP)
gen GDPperPopulation=GDP/population
gen lastGDPgrowth=L.GDP/L.L.GDP-1
gen lastGDPgrowth1=L.GDP/L.L.GDP
cap drop ten
gen ten=0
replace ten=1 if (year>=2011) & !missing(ten)
gen CO2percitylight=CO2/citylight
gen fixassetperGDP=fixasset/GDP
gen fixassetperpopulation=fixasset/population
*gen citylightbefore=ln(L.citylight)
replace citylight=citylight/10^4
*rename 第二产业占GDP比重 industry2GDP 
winsor2 pressure, cut (1 99)
winsor2 pressure2, cut (1 99)


***********************
*描述性统计
***********************
replace  expectGDPgrowth=expectGDPgrowth*100  //目标变成了100
replace  CO2perGDP=CO2perGDP*100  //每万元gdp的co2排放量


rename CO2perGDP Carbon
rename expectGDPgrowth TargetGDP
rename lastGDPgrowth GDPgrowth 
rename industry2GDP Industrialratio
rename FDIGDP FDIratio
rename GDPperPopulation GDPper

global controls "GDPgrowth  Industrialratio FDIratio  GDPper "


logout, save(描述性统计_score1) excel replace:       ///  
tabstat Carbon TargetGDP $controls    citylight pressure,  ///
stats(mean sd min p25 p50 p75 max count) c(s) 


***********************
*基准结果
***********************
*基准

reg Carbon TargetGDP
est store m_0

xtreg Carbon TargetGDP,fe
est store m_1

*加入控制变量
xtreg Carbon TargetGDP $controls ,fe
est store m_2


*cluster省份
xtreg Carbon TargetGDP $controls ,fe cluster(province)
est store m_3


*去掉北京上海广州深圳重庆一线城市
preserve
drop if city=="北京市" | city=="上海市"  | city=="广州市"  | city=="深圳市"  
xtreg Carbon TargetGDP $controls ,fe
est store  m_4
restore

*去掉北上广深一线城市   
preserve
xtbalance, range(2002 2017)  miss(Carbon TargetGDP $controls)
xtreg Carbon TargetGDP $controls ,fe
est store  m_5
restore


gen lnco2=ln(Carbon+1)
xtreg lnco2 TargetGDP $controls ,fe
est store m_6


     esttab m_0 m_1 m_2 m_3 m_4 m_5 m_6    using Table1.csv, replace    ///
     scalar( r2_w  N  F ) nogaps compress   ///
         star(* 0.1 ** 0.05 *** 0.01)  ///
         b(%6.3f) se(%6.3f)      ///
         mtitles( "1"  "2" "3" "4" "5" "6" "7")   

****************************************
*对不可观测因素的选择性偏误
****************************************
*合并市级数据

sort city year
cap drop _merge
sort city year
merge city year using 市级控制变量.dta


replace 规模以上工业企业数=ln(规模以上工业企业数)


* 加入其他可能的控制变量
**经济的 规模以上工业企业数 城乡居民储蓄年末余额

global factor1  固定资产投资总额 规模以上工业企业数
foreach v of varlist 固定资产投资总额 规模以上工业企业数{
    replace `v'=ln(`v')
}
xtreg Carbon TargetGDP $controls  $factor1,fe  


preserve
   global xvars "$controls $factor1" 
   global treat "TargetGDP"
   global depvar "Carbon"  
   xtreg $depvar $treat $xvars,fe
   scalar coef_riots = _b[$treat]
   xtreg $depvar $xvars,fe
   predict xgamma
   scalar VE=e(rmse)^2
   qui xtreg $treat xgamma,fe
   scalar coef1=_b[xgamma]
   qui xtreg $treat $xvars,fe
   predict treatres, res
   qui sum treatres
   scalar VE1=r(Var)
   scalar bias1=coef1*VE/VE1
   scalar sratio1 = coef_riots/bias
   di bias, coef_riots1, sratio1
restore



sum 人口密度 人口自然增长率 第二产业从业人员比重
**劳动力人口的 人口密度 人口自然增长率 第二产业从业人员比重
global  factor2 "人口密度 人口自然增长率 第二产业从业人员比重" 

xtreg Carbon TargetGDP $controls $factor1 $factor2,fe  

preserve
   global xvars "$controls $factor1 $factor2" 
   global treat "TargetGDP"
   global depvar "Carbon"  
   xtreg $depvar $treat $xvars,fe
   scalar coef_riots = _b[$treat]
   xtreg $depvar $xvars,fe
   predict xgamma
   scalar VE=e(rmse)^2
   qui xtreg $treat xgamma,fe
   scalar coef1=_b[xgamma]
   qui xtreg $treat $xvars,fe
   predict treatres, res
   qui sum treatres
   scalar VE1=r(Var)
   scalar bias2=coef1*VE/VE1
   scalar sratio2 = coef_riots/bias
   di bias2,  sratio2
restore




**科技的 RD人员 发明专利授权数 每万人在校大学生数 教育支出
global  factor3 " 普通高等学校学校数  教育支出" 
foreach v of varlist RD人员 发明专利授权数  教育支出{
    replace `v'=ln(`v')
}

xtreg Carbon TargetGDP $controls  $factor1 $factor2  $factor3,fe  

preserve
   global xvars "$controls $factor1 $factor2" 
   global treat "TargetGDP"
   global depvar "Carbon"  
   xtreg $depvar $treat $xvars,fe
   scalar coef_riots = _b[$treat]
   xtreg $depvar $xvars,fe
   predict xgamma
   scalar VE=e(rmse)^2
   qui xtreg $treat xgamma,fe
   scalar coef1=_b[xgamma]
   qui xtreg $treat $xvars,fe
   predict treatres, res
   qui sum treatres
   scalar VE1=r(Var)
   scalar bias3=coef1*VE/VE1
   scalar sratio3 = coef_riots/bias
   di bias3, coef_riots, sratio3
restore







***********************
*工具变量回归 借鉴nunn的方法，通过构造所在省份地级市数量( 与个体变化有关) 
*与未来两期国家经济增长目标的均值( 与时间有关) 的交互项，作
*为地级市经济增长目标约束的工具变量。
***********************
*考虑内生性 工具变量1 省的gdp增长目标
*生成每个省份中市的数量
cap drop _m
sort province
merge province using cityquantity

egen cityid=group(city)
tsset cityid year
cap drop iv
*和国家目标匹配

sort year
cap drop _m
merge  year using 国家目标
tsset cityid year
gen iv=cityquantity*F2.gdpe_n

xtivreg2 Carbon $controls (TargetGDP = provinceexpectGDPgrowth),fe
est store  m_1
xtivreg2 Carbon $controls (TargetGDP = iv),fe
est store  m_2

xtivreg2 Carbon $controls (TargetGDP = provinceexpectGDPgrowth iv),fe
est store  m_3
xtoverid

     esttab m_1 m_2 m_3  using Table1.csv, replace    ///
     scalar( r2_a  N  widstat idstat ) nogaps compress   ///
         star(* 0.1 ** 0.05 *** 0.01)  ///
         b(%6.3f) se(%6.3f)      ///
         mtitles( "IV:Provincetarget"  "IV:Num.city" "IV:Provincetarget,Num.city") 



*xtreg Carbon pressure2_w,fe

*xtreg Carbon pressure_w industry2GDP industry3GDP  FDIGDP GDPperPopulation lastGDPgrowth fixasset ,fe
*est store m_2

/*
gen crossest=TargetGDP*east
xtreg Carbon  TargetGDP crossest   ,fe
est store m_1
xtreg Carbon  TargetGDP crossest $controls ,fe
est store m_2
*/

***********************
*十三五规划
***********************
cap drop crossten
gen crossten=TargetGDP*ten 
xtreg Carbon  TargetGDP crossten   ,fe
est store m_1
xtreg Carbon  TargetGDP crossten    $controls, fe
est store m_2

     

*/ 


*xtreg Carbon  pressure2 ,fe

*xtreg Carbon  pressure_w   $controls ,fe

*xtreg Carbon  pressure_w pressure_w*district industry2GDP industry3GDP lastrealGDP FDIGDP GDPperPopulation lastGDPgrowth fixasset ,fe
 
*xtreg Carbon  pressure_w  pressure_w*ten industry2GDP industry3GDP lastrealGDP FDIGDP GDPperPopulation lastGDPgrowth fixasset ,fe



***************************
*词频分析
***************************
cap drop _m
sort city year 
merge city year using  numwords 

gen crosswords=TargetGDP*numwords

xtreg Carbon  TargetGDP crosswords   ,fe
est store m_5

xtreg Carbon  TargetGDP crosswords   $controls ,fe
est store m_6

 esttab m_1 m_2  m_5 m_6  using Table2.csv, replace    ///
     scalar( r2_w  N  F ) nogaps compress   ///
         star(* 0.1 ** 0.05 *** 0.01)  ///
         b(%6.3f) se(%6.3f)      ///
         mtitles( "1"  "2" "3" "4" )  
 


***********************
*市场机制
***********************

gen emissiontradingcity2=.
 replace emissiontradingcity2=0  if (province!="广东省" & province!="湖北省" & province!="北京" & city!="深圳市" & province!="上海" & province!="天津" & province!="重庆市" | year<2011) 
 replace emissiontradingcity2=1  if (province=="广东省" | province=="湖北省" | province=="北京" |city== "深圳市" |province== "上海" |province== "天津" | province=="重庆市") & year>=2011
 
gen crosstrade=TargetGDP*emissiontradingcity2

xtreg Carbon  TargetGDP crosstrade   ,fe
est store m_31

xtreg Carbon  TargetGDP crosstrade GDPgrowth  Industrialratio FDIratio  ,fe
est store m_41


*与碳排放交易权的数量交乘
gen city2=substr(city,1,6)
sort city2 year
cap drop _m
merge city2 year using citytrade
gen carbontrade=年碳交易总量
drop 年碳交易总量


cap drop province2
gen province2=substr(province,1,6)
sort province2 year
cap drop _m
merge province2 year using provincetrade 
replace carbontrade=年碳交易总量 if carbontrade==. &  年碳交易总量!=.

replace carbontrade=0 if carbontrade==.
gen lncarbontrade=ln(1+carbontrade)

gen crosstrade2=TargetGDP*lncarbontrade

xtreg Carbon  TargetGDP crosstrade2   ,fe
est store m_32

xtreg Carbon  TargetGDP crosstrade2  GDPgrowth  Industrialratio FDIratio  ,fe
est store m_42

 esttab m_31 m_41 m_32 m_42  using Table5.csv, replace    ///
     scalar( r2_w  N  F ) nogaps compress   ///
         star(* 0.1 ** 0.05 *** 0.01)  ///
         b(%6.3f) se(%6.3f)      ///
         mtitles( "1"  "2" "3" "4" )  




gen post=0
replace post=1 if year>=2014 & city== "深圳市"
replace post=1 if year>=2013 & province== "上海"
replace post=1 if year>=2013 & province== "北京"
replace post=1 if year>=2013 & province== "天津"
replace post=1 if year>=2013 & province== "广东省"
replace post=1 if year>=2014 & province=="重庆市"
replace post=1 if year>=2014 & province=="湖北省"
replace post=1 if year>=2016 & province=="福建省"



gen did=post*TargetGDP 
xtreg Carbon  TargetGDP did  GDPgrowth  Industrialratio FDIratio  ,fe


************************
*利用灯光做稳健性检验
************************
*利用灯光作为gdp的替代变量
replace CO2percitylight=CO2percitylight*10^4
xtreg CO2percitylight TargetGDP $controls  ,fe
est store m_1

xtreg CO2percitylight TargetGDP crossten $controls, fe
est store m_2

xtreg CO2percitylight  TargetGDP crosstrade $controls ,fe
est store m_3


 xtreg CO2percitylight  TargetGDP crosswords $controls    ,fe
est store m_4

 esttab m_1 m_2 m_3 m_4  using Table3.csv, replace    ///
     scalar( r2_w  N  F ) nogaps compress   ///
         star(* 0.1 ** 0.05 *** 0.01)  ///
         b(%6.3f) se(%6.3f)      ///
         mtitles( "1"  "2" "3" "4" )  

*********************************
*分位数回归
*********************************
cap drop cityid
egen cityid=group(city)
mat result1 = J(1,1,.) //生成空矩阵result1,2行1列
tempname a b b1
forvalues quantile = 0.1(0.1)0.9 {   //从0.1-0.9没隔0.1开始循环 
    xtqreg  Carbon  TargetGDP $controls , i(cityid) q(`quantile')   //面板分位数回归
    matrix `a' = r(table)   //对应各个分位点估计的`x1'系数
    matrix `b' = `a'[1,1] \ `a'[2,1] //仅取出面板分位数回归`a'矩阵的前两行构成`b'矩阵,第一行是系数。第二行是标准误，第四行是p值
    mat list `b'   //列出`b'矩阵  
    mat `b1' = `b''  
   matrix colnames `b1' = b_q p_q   //`b1'为`b'矩阵的转置
   mat list `b1'   //列出`b1'矩阵
 matrix result1 = (result1, `b1')    //`b'数据导入列出空矩阵result1
}

mat list result1
matselrc result1 cresult, c(2/19) //19列 = 2*分位个数+1列无用项（第一列），保留有用数据列
mat list cresult

svmat cresult, names(reg)   //矩阵数据转换为变量
keep reg*

format reg* %4.3f   //设定格式为保留三位小数
drop if reg1==. //删除缺失值，行数为8,8个解释变量：变量为18，对应不同分位点的系数和标准误；

forv i = 1/9 {  /*9：分位个数*/
local k1 = 2*`i' - 1   //奇数列为系数
local k2 = 2*`i'       //偶数列为标准误
gen up_`i' = reg`k1' + 1.96*reg`k2'   //上界
gen low_`i' = reg`k1' - 1.96*reg`k2'  //下界
}

preserve 
keep reg1 reg3 reg5 reg7 reg9 reg11 reg13 reg15 reg17   //保留系数列
rename (reg1 reg3 reg5 reg7 reg9 reg11 reg13 reg15 reg17) (reg1 reg2 reg3 reg4 reg5 reg6 reg7 reg8 reg9)   //重命名
xpose,clear  //转置数据
renvars v*,  prefix(coef_)   //重命名加前缀coef_
save coef.dta,replace  //储存系数数据
restore   //preserve+restore运行不保存

keep up* low*   //保留上下界数据
xpose,clear   //转置数据
preserve
keep if mod(_n,2)==1   //保留奇数行作为上界数据
renvars v*,  prefix(up_)     //重命名加前缀up_
save up.dta,replace    //储存上界数据
restore

preserve   
keep if mod(_n,2)==0   //保留偶数行作为下界数据
renvars v*,  prefix(low_)   //重命名加前缀low_
save low.dta,replace   //储存下界数据
restore 

use up.dta,clear    //打开上界数据
merge 1:1 _n using low.dta   //一对一合并数据
keep if _merge==3   //合并匹配成功为_merge==3，保留合并成功行
drop _merge  //删除合并信息
merge 1:1 _n using coef.dta   //一对一合并数据
keep if _merge==3
drop _merge  //删除合并信息

matrix input myrvec = (0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90)
mat list myrvec
mat qq = myrvec'  //qq矩阵为myrvec的转置
mat list qq
svmat qq, names(qqp)  //矩阵数据转换为变量

forv i = 1/1{
twoway (line up_v`i' qqp,lc(black*1.1) lpattern(dash) xlabel(0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90)) ///
(line low_v`i' qqp,lc(black*1.1) lpattern(dash) xlabel(0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90))  ///
(line coef_v`i' qqp, lc(black*1.5)  lw(thick)) ,ylabel(-0.16(0.04)0.23) yline(0, lcolor(black*2) lpattern(longdash)) xtitle("Quantile") ytitle("Effect of growth target") plotregion(margin(zero)) graphregion(color(white)) legend(off)
}




**********************************
*机制研究
**********************************
/*
clear
use "E:\论文相关\2000-2013newaa（行业小类+两高一剩）.dta"
rename 年份 year
rename 市 city
sort city year
tab 两高一剩
bysort city year 两高一剩: egen Numnewfirm=total( 新开工的企业数量 )
duplicates drop city year 两高一剩,force
keep city year 两高一剩 Numnewfirm
reshape wide Numnewfirm ,i( city year ) j( 两高一剩 )
replace Numnewfirm0=0 if Numnewfirm0==.
replace Numnewfirm1=0 if Numnewfirm1==.
sort city year
save "E:\论文相关\newfirmnum.dta",replace
*/

*被解释变量 污染企业开工率

sort city year
cap drop _m
merge  city year using newfirmnum
xtreg Numnewfirm1 TargetGDP  ,fe
est store m_1
*xtreg Numnewfirm1 TargetGDP $controls ,fe



cap drop meanNumnewfirm1
bysort 省代码 year: egen meanNumnewfirm1=mean(Numnewfirm1)
cap drop Numnewfirmraito
gen Numnewfirmraito=Numnewfirm1/meanNumnewfirm1
*gen Numnewfirmraito2=Numnewfirm1/Numnewfirm0

xtreg Numnewfirmraito TargetGDP    ,fe
est store m_2
 
*被解释变量  环境治理投资比率
cap drop _m
sort province year
merge province year using investGDP 
cap drop _m
sort province year
merge province year using invest 

gen envinvest=ln(环境污染治理投资总额亿元+1)

xtreg envinvest TargetGDP  ,fe
est store m_3

xtreg 环境污染治理投资占GDP比重 TargetGDP  ,fe
est store m_4


 esttab m_1 m_2 m_3 m_4  using mechanism.csv, replace    ///
     scalar( r2_w  N  F ) nogaps compress   ///
         star(* 0.1 ** 0.05 *** 0.01)  ///
         b(%6.3f) se(%6.3f)      ///
         mtitles( "1"  "2" "3" "4" ) 
